function [steady_state_dist,c_vec] = steadystate_distribution(bParams,oParams,f,C,lambda_mat,equity_mat,optimal_coupon,n_e,product_distribution_entrants,coupon_mat)
% Calculate the steady state distribution for firms that are in default

% Variables:
% bParams oParams   see readme
% f                 rate of creative destruction
% C                 vector with available coupons
% lambda_mat        R&D intensity 
% equity_mat        equity value
% optimal_coupon    optimal coupon
% n_e               maximum number of products entrants
% product_distribution_entrants     product distribution entrants
% coupon_mat        coupon choice

% Set parameters
	pbar = oParams.pbar;
	n = oParams.n;
    theta = bParams.theta;
    
% Find coupon elements used in model
    c_vec = unique(optimal_coupon);
    for i = 1:length(C)
        c_vec = unique([c_vec,(coupon_mat(i,:)>C(i)) .* coupon_mat(i,:)]);
    end
    c_vec = sort(c_vec);
    if sum(C==0)==0 && sum(c_vec==0)==1 % Remove zero element if its not in C
        c_vec = c_vec(2:end);
    end
    n_c = length(c_vec);

% Default thresholds
    p_D = zeros(n_c,1);
    for i = 1:n_c
        p_D(i) = sum(equity_mat((C==c_vec(i)),:)==0) - 1;
    end
    
% Set variables
    Transition_matrix =  eye((pbar+1)*n_c);
    Inflow_vector = zeros((pbar+1)*n_c,1);
    
% Loop over coupon
    for i = 1:n_c
    
% Lambda
        lambda_dum = lambda_mat(C==c_vec(i),:);

% Non-default states
        for j = p_D(i)+1:pbar
            
% Outflow (Creative destruction and innovation)
            Transition_matrix((pbar+1)*(i-1)+j+1,(pbar+1)*(i-1)+j+1) = -f * j - lambda_dum(j+1) * (1-binocdf(0,n,theta));
            if j == pbar
                Transition_matrix((pbar+1)*(i-1)+j+1,(pbar+1)*(i-1)+j+1) = -f * j;
            end
    
% Inflow
% Creative destruction (into i)
            if j < pbar
                Transition_matrix((pbar+1)*(i-1)+j+1,(pbar+1)*(i-1)+j+2) = f * (j+1);
            end
    
% Innovations (out of i)
            if j < pbar
                for l = 1:n
                    Transition_matrix((pbar+1)*(i-1)+min(j+1+l,pbar+1),(pbar+1)*(i-1)+j+1) = Transition_matrix((pbar+1)*(i-1)+min(j+1+l,pbar+1),(pbar+1)*(i-1)+j+1) + lambda_dum(j+1) * binopdf(l,n,theta);
                end
            end
        end
    
% Default (and transition to another c)
        if p_D(i) == 0
            Transition_matrix((pbar+1)*(i-1)+1,(pbar+1)*(i-1)+1) = 1;
        else
            c_dum = optimal_coupon(p_D(i)+1);
            Transition_matrix((pbar+1)*(find(c_vec==c_dum)-1)+p_D(i)+1,(pbar+1)*(i-1)+p_D(i)+2) = f * (p_D(i)+1);
        end
    
    end

% Inflow vector
    for i = 1:n_e
        Inflow_vector((pbar+1)*(find(c_vec==optimal_coupon(i+1))-1)+i+1) = -1 * product_distribution_entrants(i);
    end

% Loop to calculate steady state distribution and inflowsize
    s  =  1;
    steady_state_dist = 0;
    while abs(sum(steady_state_dist)-1) > 10^-5
    
% Update inflow vector
        Inflow_vector = s * Inflow_vector;
    
% Calculate steady state distribution
        steady_state_dist = Transition_matrix\Inflow_vector;

% Update s
        s = s/sum(steady_state_dist);
    end

end
